En la unidad anterior se describio como primer paso de los experimentos de simulación, la selección de una fuente de números aleatorios y la obtención de los mismos. Por lo tanto, disponer de un buen generador de números aleatorios es fundamental en Simulación. De hecho, constituye también una parte fundamental en otras áreas, como la Informática (algoritmos aleatorizados, verificación de algoritmos, complejidad de algoritmos, criptografía, ...), la Estadística (métodos de muestreo y remuestreo, contrastes Montecarlo, Inferencia Bayesiana, ...), y, en general, en cualquier problema de cálculo científico que, de manera directa o indirecta, incluya algún elemento probabilístico incluyendo los juegos de video.
Actualmente existen multiples generadores de números aleatorios en diferentes entornos y compiladores lo cual supondría para un usuario de la Simulación que no es necesario su estudio. Sin embargo, estudios sobre algunos generadores comerciales sugieren que debemos actuar con cuidado con el uso de ellos. Incluso, el uso progresivo de modelos de simulación cada vez más detallados exige generadores de números aleatorios de mayor calidad.
Normalmente, pasaremos un conjunto de pruebas cuando introduzcamos un nuevo generador, o cuando deseamos analizar si un generador del que desconozcamos su algoritmo subyacente es suficientemente bueno para nuestros propósitos.
Consideramos el caso en que $F_0$ es continua. La función de distribución empírica de una muestra $X_1, X_2,...,X_n$ se define como:
$$F_n(x)=\frac{\#\{X_{i}\leq x\}}{n}$$Bajo la hipótesis nula $H_0:F_X(x)=F_0(x)$, esperamos que $F_n$ se aproxime a $F_0$. Definimos el estadístico bilateral de Kolmogorov-Smirnov
$$D_{n}=\underset{x\in\mathbb{R}}{\sup}\left|F_{n}(x)-F_{0}(x)\right|$$La distribución exacta de $D_n$ está tabulada para valores seleccionados de $n\leq40$ y del nivel de significación $\alpha$. Para muestras grandes, se utiliza la distribución asintótica de $D_n$, que viene dada, para todo $z\leq0$, por
$$\underset{n\rightarrow\infty}{\lim}P(\sqrt{n}D_{n}\leq z)=L(z)=1-2\overset{\infty}{\underset{i=1}{\sum}}(-1)^{i-1}e^{-2i^{2}z^{2}}$$$L(z)$ está tabulada y se comprueba que la aproximación es suficientemente buena para $n\geq35$. Intuitivamente, esperamos que $D_n$ sea pequeño cuando la hipótesis nula es cierta. En nuestro caso particula de aleatoriedad, si $X_{(1)}<X_{(2)}<...<X_{(n)}$ designa al estadístico de orden, $F_0(x_{(i)})=x_{(i)}$, y como $F_n(x_{(i)})=i/n$, resulta:
$$D_{n}=\underset{i\leq i\leq n}{\max}\left\{ \max\left[\left|\frac{i}{n}-x_{(i)}\right|,\left|x_{(i)}-\frac{i-1}{n}\right|\right]\right\}$$
In [12]:
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
N = 100
In [13]:
# Teoretical CDF for Uniform distribution F(x)
u = [1.0] * N
#print(u)
# Gerate random numbers with Uniform distribution f(x)
#np.random.seed(123456789)
x = np.random.uniform(0,1,N)
# Plot empirical distribution f(x)
count, bins, ignored = plt.hist(x, 25, normed=True)
plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
plt.show()
In [14]:
x1 = np.arange(1/N,1+1/N,1/N)
y1 = np.cumsum(np.sort(u)/np.max(np.cumsum(u)))
x2 = np.sort(x)
y2 = np.cumsum(np.sort(x)/np.max(np.cumsum(x)))
De=np.absolute(y2-y1)
#print(D)
print("De = ",np.max(De))
# Plotting
plt.plot(x1,y1,label='Uniform')
plt.step(x2,y2,label='Random')
x3 = [x1[De.argmax()],x1[De.argmax()]]
y3 = [np.min([y1[De.argmax()],y2[De.argmax()]]),np.max([y1[De.argmax()],y2[De.argmax()]])]
plt.step(x3,y3,label='D')
#plt.plot(x1,De,label='De')
plt.legend()
plt.xlabel('X')
plt.ylabel('Cumulative Probability P(X)')
#plt.title('arbitrary')
#plt.xscale('log')
#plt.yscale('log')
#plt.xlim(0, np.max([np.max(whitney),np.max(redwell)])+0.05)
plt.ylim(0, 1.05)
plt.show()
In [16]:
# scipy function for K-S test
from scipy import stats
D,pvalue = stats.kstest(y2, 'uniform')
print("D = ", D)
print("p-value = ", pvalue)
La prueba $\chi^{2}$ es de bondad de ajuste. Es poco potente, por lo que permite justificar el rechazo de una hipótesis, pero proporciona escaso soporte a su aceptación. El problema de bondad de ajuste se plantea como sigue. Tenemos una muestra $X_1, X_2,...,X_n$ de una población con distribución $F_X(x)$ desconocida. Deseamos contrastar la hipótesis nula $H_0:F_X(x)=F_0(x)$, para todo $x\in\mathbb{R}$, donde $F_0(x)$ está completamente especificada, frente a la alternativa $H_1:F_X(x) \neq F_0(x)$ para algún $x$. Para realizar el contraste, partimos el soporte de $X$ en $k$ subconjuntos o clases mutuamente excluyentes. Para cada subconjunto $i$, $i=1,...,k$, calculamos el número $f_i$ de observaciones que caen en la clase y el numero esperado de observaciones $e_i$ bajo la distribución hipotética $F_0$. El estadístico de la prueba es:
$$\chi^{2}=\overset{k}{\underset{i=1}{\sum}}\frac{(f_{i}-e_{i})^{2}}{e_{i}}$$cuya distribución asintótica es una $\chi_{k-r-1}^{2}$, donde $r$ es el número de parámetros que se estiman a partir de las observaciones. La aproximación se considera aceptada si $\min\left\{ e_{i}>5\right\}$. Intuitivamente, el estadístico tiende a ser pequeño bajo la hipótesis nula, y grande cuando ésta es falsa, con lo qu ese deduce inmediatamente un contraste de hipótesis. La potencia del contraste crece con el número $k$ de subconjuntos considerados.
En nuestro problema de aleatoriedad, $F_0$ es la función de distribución de la uniforme en [0,1], con lo que $r=0$. Además, pueden cogerse $k$ subintervalos de [0,1] de igual longitud, con lo que $e_i=n/k$.
Dada la sucesión de observaciones $X_1, X_2,...,X_n$, construimos la sucesión de simbolos binarios definida mediante 1 si $X_i<X_{i+1}$, 0 si $X_i>X_{i+1}$. Definimos racha creciente (decreciente) de longitud $l$ a un grupo seguido de $l$ números 1 (ó 0). Contabilizamos el número de rachas. Sabemos que su distribución asintótica, bajo la hipótesis nula de aleatoriedad, es:
$$N\left(\frac{2n-1}{3},\frac{16n-29}{90}\right)$$Intuitivamente, rechazamos la aleatoriedad con un número muy pequeño o muy grande de rachas. De ahí se obtiene inmediatamente la prueba.
Otro procedimiento para definir rachas se tiene del recuento de observaciones que se sitúan a un mismo lado de la mediana. La distribución asintótica del número de rachas, bajo la hipótesis de aleatoriedad, es:
$$N\left(1+\frac{n}{2},\frac{n}{2}\right)$$de donde se sigue, inmediatamente una prueba.
Separamos las observaciones en $k$-uplas:
$$(U_1,U_2,...,U_k),...,(U_{ik+1},U_{ik+2},...,U_{(i+1)k}),...$$La $k$-upla general se escribe:
$$(U_{ik+j})_{j=1}^k$$La ordenamos crecientemente y consideramos la ordenación correspondiente de los subíndices $j$. Bajo la hipótesis de que la probabilidad de que dos números sean iguales es nula, hay $k!$ ordenaciones posibles. Bajo la hipótesis de independencia, todas las permutaciones son equiprobables, con probabilidad $1/k!$. Entonces es inmediato aplicar una prueba $\chi^2$ con $k!$ clases, distribución asintótica $\chi_{k!-1}^2$, frecuencias esperadas $r/k!$, onde $r$ es el número de $k$-uplas y frecuencias observadas el número de veces que aparece cada ordenación.
Fijamos dos valores $\alpha$ y $\beta$ con $0<\alpha<\beta<1$. La sucesión presenta un hueco de longitud $m$ si $U_j, U_{j+m}\in\left[\alpha,\beta\right]$ pero $U_{j+1},...,U_{j+m+1}\notin\left[\alpha,\beta\right]$. Bajo la hipótesis de aleatoriedad de la serie, la longitud $m$ de los huecos sigue una distribución geométrica de parámetro $P(\alpha\leq X_{j}\leq\beta)=\beta-\alpha=\rho$, es decir:
$$P(hueco\: longitud\: m)=\rho(1-\rho)^{m-1}$$.
La hipótesis de aleatoriedad implica independencia de las longitudes de los huecos y podemos aplicar una prueba $\chi^2$ basado en las comparacioens de los números observados y esperados de huecos de longitud $m$.
In [ ]: